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The N-body problem in which the force on the Ith body is 
Fj = k ^2 j mjmj(xj — x/)/r 3 where r is the time dependent radius of the whole 
system is solved exactly for all initial conditions. Here, writing x for the position 
of the barycentre, r 2 = M~ l J2i m l{' x -l ~ x) 2 . 

Also solved are all N-body problems for which the potential energy is of the 
form V = V(r) for any function V where r is given above. The first problem is the 
special case V oc r _1 . In all cases r vibrates for ever, just like the radius to a body 
which has angular momentum about a centre of attraction due to the potential 
V(r). Violent Relaxation does not occur in these systems. When r vibrates as 
in a highly eccentric orbit the exact solution should prove to be a useful test of 
N-body codes but a severer test follows. 

We show further that when Fj is supplemented by an inverse cubic repulsion 
between each pair of bodies —k'J2j^i m i m j( x J ~ *-i)/(xj — x i) A there is still no 
damping of the fundamental breathing oscillation which performs its non-linear 
periodic motion for ever. This is a special case of continually breathing systems 
with potentials Vo(r) + r~ % Vi ((xj — xj)/r). 

A remarkable new statistical equilibrium is found for systems whose size is 
continually changing! 

Keywords: N-body systems; Newton; orbits; dynamical systems; 
statistical mechanics; pulsating equilibria 



1. Introduction 

In 1995 one of us (Lynden-Bell, R.M. 1995, 1996) showed that the classical 
statistical mechanics of systems of N bodies interacting through mutual potential 
energies of the form 

V = NF(J2<l 2 i/N) (1.1) 

could be exactly calculated for any chosen function F. By choice of a suitable 
function F, she found a case with a simple phase transition which was calculated 
for any N. 

When the particles are of equal mass and when q/ measures the displacement 

of the I th particle from the barycentre x, J^Q/ /N may be rewritten as a mutual 

I 
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potential energy 

iV-^qf = 7V- 1 ^(x J -x) 2 = £ E( X /- X ^V 2 • (1-2) 
/ I < J 

More than 300 years previously Newton showed that the N-body problem in 
which the force on body I due to body J was 

km I m J (y. J - xj) 

could be exactly solved with all the interactions present (Newton (1687), see also 
Chandrasekhar (1995)). Newton's case corresponds to a potential energy 

V = h k J2J2 rnjmj^j - xj) 2 . (1.3) 
I < J 

When the masses are equal the expression (1.3) is the linear- F case of expression 
(1.1) by virtue of (1.2). As the statistical mechanics associated with expression 
(1.1) is still beautiful even when F is any non-linear function, we are led to wonder 
whether Newton's exact solution to the N-body problem could be generalised to 
cover potential energies of the form V = V(r). Here V is any function and r is 
the mass weighted root mean square size of the system 

r 2 = M- 2 E mimj^! - xj) 2 = M'^mi^i - x) 2 (1.4) 
I < J I 

where 

M = E m -f an d x = M~ 1 E m -f x / • (1-5) 

Following this idea we found a beautiful dynamics in 3iV dimensional space 
which precisely mimics normal orbit theory in 3 dimensions. 

2. Solutions to the N-body Problems 

The equations of motion are 

m/ x 7 = -dF/dx/ . (2.1) 

Provided V is a mutual potential energy, i.e., a function of the x/ — x j then 

XdF/dx/ = so X = and 
I 

x = x + ut (2.2) 
Writing q/ = x/ — x we have r 2 = M -1 J]m/q 2 so in centre of mass coordinates 

the equations of motion for the q/ are 

dV dV mi qj 

m ™ = ~d^ = Tm7 • (2 - 3) 

We now invent a 3iV dimensional vector r of length r, whose components are 

M qi ' VM q2 VM q/ ■■■)■ {2A) 
The first 3 components of r involve the position of particle 1 and the next 3 
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particle 2, etc. Multiplying (2.3) by ^M/mj we see that those equations for all 
I can be rewritten using unit vector r 

dV r 

Mr = -—- = -V'r. (2.5) 
or r 

For insight please consult here the note added in proof. 

Let a and (5 be indices that run from 1 to 3N. Then taking the antisymmetric 
product of (2.5) with r we find 

M ~<ti ~ T ^ aS> = ° 

so the hyper-angular-momenta (per unit mass) are all conserved 

r a r/3 - rpr a = L afi = const = —Lp a . (2.6) 

These quantities are not just normal angular momenta. They involve among 
others the x component of particle 1 multiplied by the rate of change of the x 
component of particle 2, i.e., x±X2 and also such terms as xii/2, etc. Only when 
both q and (5 are in the range 1, 2, 3 does L a p reduce to components of the 
normal angular momentum of particle 1 (in that case). The trace of the square 
of the antisymmetric tensor L a p is 



LapLap — 2 



r l v z - (r • if = 2L 2 . (2.7) 



Notice that if we were in 3 dimensions L? would be just (r x r) 2 so it would 
then be the square of the specific angular momentum. We shall refer to ML as 
the total hyper-angular-momentum and ML a p as the hyper-angular-momentum. 
Taking the inner product of equation (2.5) with r and integrating, we obtain the 
energy equation 

\Mi 2 + V{r) = E . (2.8) 

Now 

r= ^(r-r) 1/2 = r • r/r (2.9) 

so from (2.7) 

r 2 = LV 2 + r 2 (2.10) 

and (2.8) may now be rewritten 

\M(r 2 + L 2 r- 2 ) + V = E . (2.11) 

Differentiating (2.11) and dividing by f we get the equation of motion for scalar 
r 

M(r - L 2 r~ 3 ) = -V'(r) . (2.12) 

This we recognise as the equation of motion of a particle of mass M and angular 
momentum ML in the central potential V(r). Equation (2.11) is of course its 
energy equation and it may be integrated in the form 

f r dr 
t ~ t ° = J y/2M-L(E-V(r))-L*r-* ' ^ 
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which gives the standard relation for the periodic]]] function r(t). Equation (2.13) 
integrates nicely for the standard cases V oc 1/r, V oc r 2 and the isochrone 
V oc (6 + s) _1 where s 2 = r 2 + b 2 . These potentials may be supplemented by 
an extra term in r~ 2 giving a r~ 3 supplementary force but they comprise the 
only central orbits that can be solved involving nothing more complicated than 
trigonometric functions (Eggen et al. 1962, Evans et al. 1990). There are many 
others that involve elliptic functions, etc., but for them the expression (2.13) is 
as simple as the function. We shall return to (2.13) later but for the present we 
merely note that it determines the periodicf function r(t). 

By analogy with normal central orbits the form of equation (2.12) strongly 
suggests that we invent an angle (j) such that 

r 2 4> = L . (2.14) 

Then from (2.11) we eliminate t and integrate to find 

L dr 



" 2 ^J2M~ 1 [E + V{r)\ — L' 2 



r 



(2.15) 



just as in ordinary orbit theory. We define 4> by (2.15) and so there is no constant 
of integration. It gives <fi as a function of r, or r as a periodic function of (j). Thus 
we have determined scalar r but we need to find all components of the 3iV vector 
r to determine the motion completely. Although V'{r)/r may now be regarded 
as a known function of t nevertheless the form of (2.5) does not at first look 
very encouraging, except for Newton's case for which V = ^kM 2 r 2 so V'(r)/r is 
constant. However, a wonderful simplicity will emerge for the general case too. 
For a moment we discuss Newton's case. Equations limited to this case will have 
the letter N appended to the equation number. Evidently (2.5) becomes 

r = -kMr (2.16N) 

so each component of r vibrates harmonically, all with angular frequency \fk~M. 
Thus each body orbits in an ellipse centred on the barycentre and moving with it 
and the whole motion relative to the barycentre is periodic with period 2iv/^/kM. 
Newton obtained this solution by realising that the linear law of attraction, when 
suitably mass weighted, gave a net total force on each particle due to all the 
others which was directed toward the barycentre and in constant proportion to 
distance from it. 

A similar simplicity emerges in the general case when we solve for r = r/r as 
a function of <f>. To do this we write 

d/x -\ d I ' L dr A L 2 d 2 r A 

r = jt ( rr + rr ) = jt b# + rr ) = 73 w + rr ' < 2 - 17 ) 

where two terms cancelled at the last step. Now scalar r is given by (2.12) so, 
inserting (2.17) into (2.5) and multiplying by r s L~ 2 , we find the lovely result 

d 2 r/d(f) 2 + r = . (2.18) 
Thus the unit vector r vibrates harmonically in (j) with period 2?r. For Newton's 



f Periodic when the system is bound. 
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case, and for the hyper-Keplerian case V oc 1/r, the magnitude r is also periodic 
with period 2tt in (p. But more generally, as in the isochrone, r is periodic but with 
a 4> period that is incommensurable with 2tt (in general). The general solution of 
(2.18) is given by 

r = A sin + B cos (/> (2.19) 
where A and B are 3iV-vector constants of integration. But r is a unit vector so 

1 = A 2 sin 2 </> + 2A-Bsin0cos</> + B 2 cos 2 c/> = ±(A 2 + B 2 )- 

-i(A 2 -B 2 )cos2</> + A-Bsin2</> . 

Since this must be true for all (f> we find 

A 2 = B 2 = 1 and A • B = , (2.20) 

so A and B must be orthogonal unit 3iV-vectors. There is a further restriction 
on these constants; because the qj are the centre of mass coordinates they must 
obey 

= (Ev^3/-3+^ r for j = 1, 2, 3 




= ^2V™I A 3I-3+j 



and 



(2.21) 



= J2V™~I B 3I-3+j 
I 

To construct constants A and B satisfying the constraints (2.20) and (2.21) 
we proceed as follows: 

Take m\ to be the largest mass; choose B a a > 3 and A a a > 3 arbitrarily, 
but set Aj = ^/mi A 3l-3+j f° r 3 = 1: 2, 3 and a similar relation for Bj. Then 

Aj and Bj satisfy (2.21) and so will AA and /u(B — i>A) where A, n and v are 
any scalars, since (2.21) are linear. 

Now choose A = |A| _1 and set A = AA so that A is a unit vector; furthermore 
set B' = B — (B • A)A so that B' • A = 0; finally normalise by writing B = \xB! 
with jj, = p^ry . Then A and B so constructed are the general unit vectors satisfying 
the constraints (2.20) and (2.21). Thus our general solution for r is 

r = rr = r(Asin</> + Bcos0) (2.22) 

Thus relative to the centre of mass every particle's orbit becomes a centred ellipse 
after a universal scaling by r _1 . This rescaling factor is time-dependent and r 
behaves like the radius to a particle in a central orbit which is an eccentric ellipse 
in the Keplerian case. This is our main result. 

Differentiating (2.22) with respect to t one finds that (using r 2 <j) = L) 

L a p = r a rp - rpr a = (B a Ap - A a Bp)L . (2.23) 

It is of interest here to count constants of integration. A and B together have 
6iV components but (2.20) and (2.21) give nine constraints, thus there are 6A^-9 
freedoms so far. In the solutions for r(t) and r{4>) there are three further constants 
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L,E and to. Thus r{t) depends on 6./V-6 constants. To these we may add the 
constants xo and u involved in the motion of the centre of mass, and our final 
solution depends on 6N arbitrary constants. This equals the number of initial 
positions and velocities that determine the motion, which verifies that we have 
the general solution. 

The solution is only fully explicit when we can carry out the integrations (2.13) 
and (2.15) which determine r and <p. We do not need this for Newton's case 
because we already gave the solution. When V oc r _1 , the hyper-Keplerian case, 
it is actually simpler to follow Hamilton's treatment of the problem with his 
eccentricity vector later used by Runge in 1919 and Lenz in 1924 (see Hamilton 
1845, Laplace 1799, Goldstein 1976, or Chandrasekhar 1995). This was certainly 
known to Laplace and Bernoulli and can be traced to Newton. The method avoids 
the awkwardness of having to treat the elliptic, parabolic and hyperbolic cases 
separately. With V = —kM 2 /r equation (2.5) becomes 

r = -kMi/r 2 (2.24K) 

contracting with L a p we find 

d ( V R V 

j t {L a prp) = -kM [r a —rf3 - fp^-rp 
d ( r, 



= kM- . 

dt\r 

Hence 

L a/3 r p = kM(r a + e a ) (2.25K) 

where e is a 3 N- vector eccentricity. Contracting with r a /r and using the anti- 
symmetry of L a p we find, 

£/r = l + e-r, (2.26K) 

where 

i = L 2 /{kM) . (2.27K) 

If e is the magnitude of e and <j> is the angle between e and f , then this equation 
becomes the equation of a conic section in the r, (j) plane, i.e., 

£/ r = l + ecos0. (2.28K) 

To demonstrate that <p is indeed the angle called by that name previously, 
we first square (2.25K) writing out L a p in full and using (2.7) on the left and 
eliminating e • f via (2.26K) on the right. This gives just as for the 3 dimensional 

case 

L 2 r 2 = k 2 M 2 (e 2 - 1) + 2L 2 Mk r - 1 , 
which we rewrite in the form of the energy equation (2.8): 

M , 2 _tM != ,tW-l) =E (229K) 

Li 



1 i> /.f,',:' 



r 



Remembering (2.10) we have 

■£ 
r 



I E M\ L 2 _ k 2 M 2 



e 2 — I - — 1 



(2.30K) 
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but by (2.30K) and (2.28K) 



= es in0(r 2 </>) , (2.31K) 
kM 

r = esin0. (2.32K) 

L 

Thus r 2 (j) = L and our new cj> can only differ from our old one by a constant x- 
Such a constant is irrelevant as it merely makes a transformation 

A — » A cos x ~ B sin x 
B — ► B cos x + A sin x 

on our constants of the motion. The constraints (2.20) and (2.21) are invariant 
to such transformations. Thus our orbit in 3iV space is given by 

£ 



1 + e cos ( 



(Asin<^ + Bcos0) . (2.33K) 



This is our prime result and gives the x/ via (2.4) and (2.2). Relative to the 
centre of mass each particle's orbit lies in a plane and eliminating cf) in favour of 
coordinates x and y in that plane we find it is a conic section. When e < 1 the 
ellipse has in general neither its focus nor its centre at the centre of mass. This 
is clear from the geometrical interpretation given in the note added in proof. 

The relationships between r, <j> and t are given as usual by Kepler's equation 
which one may derive from (2.13) and from (2.28K) and (2.31K): 



dr 



d 



( i \ r3 

\ l+ecos<^ J 



dt = — = ^7 = ~< ' (2.34K) 

r M*esin0 k 2 M 2 (1 + e cos (j)) 2 v ' 

Defining a = GM/(—2e) these give, writing r = a(l — ecosn), 

t - t min = GM(-2e)~ 3/2 (77 - e sin rj) , (2.35K) 

and 

4> = 2 tan" 1 ^YTe 2 j ' (2.36K) 

where eft is measured from pericentre. Prom this one deduces cos <p = and 

£ = a(l — e 2 ). The generalisations of these formulae for the isochrone potential 
are given in the Appendix. 

3. Generalisation to Breathing Systems and their new Statistical 

Mechanics 

When the potential is of the form]]] Vo(r) + r~ 2 V2(r) equation (2.5) becomes 
Mr = -Vo(r)r + 2r" 3 y 2 r - r~ 3 dV 2 /dx , (3.1) 

f r is constrained by the barycentre constraint. It is understood that V2 takes values dependent only 
on the subspace to which f is constrained 
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where the last term in the force is not hyperradial. Indeed since V 2 is independent 
of r it is constant along each radial line so its gradient is automatically perpen- 
dicular to r. We therefore form a Virial theorem by taking the dot product of 
(3.1) with r and using 

r.r=| (r . r) - r 2 = I^)-r 2 . (3.2) 

We thus obtain writing T = Yl \Mr 2 

d 2 

\M-^{r 2 ) = 2T- rV ' + 2r~ 2 V 2 . (3.3) 
But by energy conservation 

T + V + r~ 2 V 2 = E , (3.4) 

so 

\M^(r 2 )=2E- l -i(r 2 V ). (3.5) 



Notice that V 2 has disappeared from this equation which is now an equation for 
the scalar r only. This separation of the motion of the scaling coordinate r from 
the rest of the dynamics only occurs when the potential is of the special form we 
have chosen with the 'angularly' dependent part scaling as r~ 2 . Multiplying (3.5) 
by dr 2 /dt and integrating 




2Er 2 - 2r 2 V - ML 2 , (3.6) 

where the final term is the constant of integration and C has the same dimensions 
as our former L whose part it plays. Dividing by 2r 2 we now have 

\M{r 2 + C 2 r- 2 ) + V Q {r) = E . (3.7) 

The equation of motion for r follows by differentiation 

M(f - £ 2 r~ 3 ) = -dV /dr . (3.8) 

Thus r vibrates as though it were in a central orbit of a particle of mass M and 
angular momentum ML with a central potential Vo(r). Integrating (3.7) we see 
that for bound orbits r vibrates for ever with a period 

* (3.9) 



^M-^E - V (r)) - C 2 r~ 2 ' 

In the above we have assumed that C 2 is positive and that Vo(r) is not so 
strongly attractive that it beats the centrifugal repulsion. Then for bound states 
there are normally simple zeros of the expression in the surd above. When the 
constant C 2 is negative, one may change the nomenclature to make it positive 
by adding a constant K to V 2 and subtracting K/r 2 from Vq so as to leave V 
unchanged. While this makes C 2 positive it may make Vq so attractive that r 
can reach the origin. If we again define an angle <\> such that r 2 cf) = C and plot 
an orbit for r in r, eft space, then such orbits approach the origin in finite time 
following logarithmic spirals (for small r). They re-emerge with (j) still increasing 
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but with r now increasing on logarithmic spirals. Although the time near the 
origin is finite, the increment in (j) on passing through the origin is infinite. The 
outside turning points of these orbits will be just as for others. 

The form of the potential energy in (3.4), Vq{t) + r~ 2 V2(r) is the general- 
isation to 3iV dimensions of the 3 dimensional potential Vq + r~ 2 V2(9,(j)) for 
which the radial motion separates from the transverse through the exact inte- 
gral ^mh 2 + V2 (#,</>) = const where h = r x v. Further separations occur if 
V2 = U(9) + W<p{4>)/ sin 2 9. Still further integrals can occur for special Vq such as 
the simple harmonic and the Kepler potential, see Evans, 1990. In 3./V-3 dimen- 
sions each separable potential gives an exactly soluble N-body problem. Marshall 
& Wojciechowski show that the general case occurs in hyperellipsoidal coordi- 
nates in 3iV-3 dimensions and that other cases, such as the hyperspherically 
separable one we have used, can be found as degenerate cases of hyperellipsoidal 
coordinates. 

We now apply the general theory of this section to the case mentioned in the 
abstract for which 

V o = \k E E ~ xj) 2 = \kM 2 r 2 , (3.10) 

1 < J 

and 

r~ 2 V 2 = \k' E m i m A*i ~ x </)~ 2 > (3- 11 ) 
1 < J 

with k and k! positive. Then we have the system of N bodies which attract each 
other linearly according to their separations and repel each other with an inverse 
cubic repulsion. Since the scale r of x/ — x j cancels out in V2 we see that indeed 
V2 = V^(r). Thus the general theory applies and the period of r is given by (3.9) 
which gives 

(3.12) 

This is only half the period in <f> because for a central ellipse r undergoes two 
oscillations as (j) increases by 2ir. 

Thus the radius of such a system will continue to vibrate for ever and shows 
no Violent Relaxation (Lynden-Bell, 1967). With both attractive and repulsive 
forces present, it may be possible to get a giant lattice solution but the increased 
pressures near the centre will give interesting radial distortions to any such lattice 
as in a planet. 

Even though the above hyper-radial motions are simple the f motions are not, 
may show relaxation phenomena and we find below that they may even achieve a 
slightly modified form of equilibrium even though the scaling variable r continues 
to vibrate at large amplitude! 

To see this we first return to the general case with Vq (r) and V2 (r) any functions 
and look for the equations of motion of the r. These follow from (3.1) if we 
remember that 

d , X . ~\ ..~ 1 d . c y\». 

r = — (rr + rr) = rr H -fr r) ; 

dt r dt 
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substituting this into (3.1) and eliminating f via (3.8) yields 

Putting r~ 2 d/dt = d/dr this may be rewritten as an autonomous equation for 
r(r) 

Md 2 i/dT 2 = -ML 2 i + 2V 2 i - dV 2 /dr . (3.13) 

The radial component of this equation is irrelevant; it follows from the condition 
|f | = 1 and the transverse components. We now write Q/ = ymi/M qj/r so 
that the Q/ are the components of r taken in triples. Then the centre of mass 
constraint becomes ^2j^/mi Q/ = 0. When we expressed V 2 in terms of our 
original x/ it depended only on differences so V 2 (xi + A, . . . xj + A . . . xn + A) 
was independent of A. This means that when we write V 2 as a function of our 
new coordinates Q/ it has the property that 

, 2 ( Ql + ^,... Qw + ^) 

is independent of A. Differentiating with respect of A and then setting A = 
we deduce that 

J2V™i dv 2/ d Qi = • ( 3 - 14 ) 
I 

Now consider the Lagrangian L(Q' / ,Q/) where Q' 7 = dQi/dr and 

L = £lMQ' / 2 -F 2 (Q 1 ...Q JV ) , (3.15) 

and the Qj are subject to the constraints 

Y^Q 2 i = l and ^v / ^7Q/ = °- ( 3 - 16 ) 

Consider also the variational principle in r-time 5 J L(Qj, Q/)dr = 0. 
Lagrange's equations are 

MQ" = -dV 2 /dQ I + A(t)Qj + /x(r)V^7 , (3.17) 

where A and n are Lagrange multipliers corresponding to the constraints (applied 
for all t). Multiplying (3.17) by yjmj and summing using (3.14), the second 
constraint equation (3.15) gives /x = 0. Evidently (3.17) and (3.13) are equivalent 
since r = (Qi, Q2 • • • Qat)- So (3.13) follows from the Lagrangian L. The 'energy' 
equation for the r follows from r 2 (3.4) — ^(3.6) when we use r 2 = r 2 (dr/dt) 2 + r 2 
and we get 

W{%) 2 + V 2 (i) = \MC 2 , (3.18) 
where by definition of Q/ 

In the statistical mechanics that follows we shall be concerned with the case in 
which N is large so that there are many terms in this sum. 
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(a) Statistical Equilibrium of a system in large amplitude oscillation 

In one sense it is surprising that a system may achieve a statistical equilibrium 
while its scale continues to oscillate (or in the unbound case evolve). However 
this only occurs for these special systems in which the large scale oscillations 
separate from the rest of the dynamics and only then when the potential Vi is 
sufficiently complicated to give quasi-ergodic motion in r. Once the r motion 
has been separated off, we have shown that the motion of Q/ follows from a 
Lagrangian in r-time and that Vi can be any function of the normalised relative 
coordinates, so it can achieve the necessary complications. Since the Q/ motion 
pays no attention to the r motion there is no reason why it should not achieve 
such an equilibrium. Whether such an equilibrium is a thermal equilibrium or not 
depends on semantics. The quantity shared by the different Q/ is not energy but 
hyper-angular-momentum, MC 2 , and the time is replaced by r but, as we shall 
see, the equilibrium remains at all phases of the oscillation and does not depend on 
the timescale of the relaxation as compared with the period of oscillation. In many 
respects it is a true generalisation of the usual concept of statistical equilibrium 
albeit with sufficient differences to make it interesting. We have N, Q/ subject to 
the 4 constraints 



and having equations of motion following the Lagrangian (3.15) and having 'en- 
ergy' given by (3.18). When there are very many particles present the fixed centre 
of mass constraints (3.19) are unimportant. They only give changes of order 1/N 
in the results and are in any case statistically satisfied by the equilibrium found 
without imposing them. As they significantly complicate the arguments while 
adding very little to the result we shall now solve the problem when only the 
constraint (3.20) is imposed and we shall specialise to an interaction V2 equiva- 
lent to a hard sphere gas. That is we shall take V2 to be negligible at any one 
time but nevertheless to be present to perform its role of redistributing hyper- 
angular-momentum. Then by (3.18) our 'energy' equation takes the form 



where on the right we have a conserved quantity. 

As far as the statistical mechanics is concerned the Q/ are equivalent, so we 
invent a 6 dimensional phase space Q, Q' and divide it into cells of equal volume. 
The number of ways W of putting n a distinguishable particles in the a th cell 
centred on Q a , Q' a is 




(3.19) 



£Q? = i 



(3.20) 




We maximise W subject to the constraints 



£n a = AT 
E*«Q? = £ 2 



(3.21) 
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and using Lagrange multipliers a, ^/3, ^7 for the constraints we obtain 
S\nW - J2 Sn a (a + ^Q' 2 + ±7Q 2 ) = . 
The normal use of Stirling's theorem then leads to 



n a = exp ■ 



a 



l + ^Q' a 2 + 7^Qal • (3-22) 



Returning to the constraint equations (3.21) we see that they are the same 
with Q' 2 /C 2 and Q 2 interchanged. Hence when we solve them for f3£ 2 and 7 
we inevitably find those to be equal. If we now pass to a continuous distribu- 
tion function /(Q, Q') in our six-dimensional phase space by replacing n a by 
f(Q,Q')d 3 Qd 3 Q' we find 

/(Q, Q') = ^expf-^Q' 2 + £ 2 Q 2 )] . (3.23) 

Where Q' 2 stands for any one of the equivalent Qf and Q 2 stands for one of the 
equivalent Q 2 = jjqj/r 2 . We shall now re-express the Qj in terms of the q/ and 
the q/ 



, 2 mi /q/V 2 r 4 

Ql= W{v) =M mi 



d_ /qj 

dt\ r 



2 ^-.2 / ■ \ 2 

= M mi ^"r q/ 



Inserting these expressions into (3.23) and writing H for the 'Hubble' expansion 
rate r/r, which is of course time dependent in our problem, the expression be- 
comes simplest written in terms of the peculiar velocity vj relative to the Hubble 
flow 

V/ = q/ - Hqj . (3.24) 
The distribution of vj and q/ at given r is then 



/(v/,q 7 |r) = A exp 



(3r 2 f l 2 \ (3C 2 Am/ 2 



(3.25) 



r is continually changing in time and by no means slowly, nevertheless $ and C 2 
are exactly constant however slowly the interactions have led to the equilibrium 
distribution. We notice that at each r the peculiar velocity vj with respect to the 
'Hubble' flow Hqj is Maxwellianly distributed with a temperature proportional 

to r~ 2 and the heavier particles move more slowly with respect to that flow. 

1/2 

Likewise m/ qi are Gaussianly distributed with a dispersion proportional to r. 

The value of is readily deduced from the condition ^qf = r 2 which gives 

(3 = (3iV-l)/£ 2 . By (3.23) / is unchanging so there is no entropy change. Here the 
— 1 follows from the fact that there are 3iV-l independent Qj components when 
the constraint (3.20) is imposed. Had we imposed also the bary centre constraints 
we would have found 37V-4. The sharing of C 2 without the sharing of energy has 
arisen before in Cometary Theory (see Rauch & Tremaine, 1996). 



4. Conclusions 

It is natural to ask whether the new form of equilibrium (3.25) can be gener- 
alised to the Hubble flow of relativistic cosmology and to look for connections 
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with the interesting fact that the cosmic background radiation remains in equilib- 
rium during cosmic expansion and creates no entropy while a classical gas creates 
entropy during uniform expansion due to its bulk viscosity. These are beyond the 
scope of this paper but we hope they show that exploration of somewhat eso- 
teric N-body problems can have considerable interest outside classical dynamics. 
Readers will have appreciated the beauty of these N-body problems treated ex- 
actly in 3N dimensions as well as the novelty of the unusual statistical mechanics 
of systems with isotropic expansions as seen in the distribution (3.25) where the 
velocities are those relative to the time-dependent flow as given by (3.24). 



Appendix A. 

Motion in Henon's (1959) isochrone potential V = —GM 2 /(b + s) where s 2 = 
r 2 + b 2 . We write e = E/M for the specific energy. We need to evaluate the 
integral (2.13) for t and (2.15) for <f>. If we write 

a = GM/(-2e) , (Al) 

and put 

(l-e 2 ) = L 2 /[-2e(a-b) 2 ] , (A2) 
then the surd in (2.13) xr is 

S = [2e{s 2 -b 2 ) + 2GM{s-b)- L 2 ] 1 ' 2 = (-2e) 1 / 2 [(a- 6) 2 e 2 - (s-a) 2 } 1 / 2 ; (A3) 
so we write 

(s — a) = — (a — b)e cos n , (A4) 

and (2.13) becomes 

t-to = J (s/S)ds = a(-2e)" 1/2 (?7 - (1 - 6/a)esin??) . (A5) 

Putting 

k = (-2e) 1/2 /a = (GM/a 3 ) 1/2 = {-2ef' 2 /GM , (A6) 
then 2-k/k is the radial period as for Kepler's case and 

n(t — to) = r] — (1 — 6/a)esin?7 , (A7) 

which we recognise as the generalisation of Kepler's equation (b = 0). Both in 
Kepler's case and for the isochrone, direct evaluation of (2.15) in terms of n gives 
untidy formulae (cf. (A17)). In the Keplerian case one uses \jr as the integration 
variable to get £/r = 1 + ecoscj). For the isochrone one has from (2.15) 

and we need to use ^ and instead of l/r. We therefore use two new angles 
X and x+ an d rewrite (A4) in terms of a p the pericentric value of s — b. 

/ 2 , 7 2 , u a p n \ Op(l + e) , (a p + 26)(l + e + ) 

y'r 2 + b 2 -b = s-b = — — (1 — ecosr/) = — — — = -26+—^ ^ 

1 — e 1 + ecosx (l + e + cosx+) 

(A9) 
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where 

a p = ( a -6)(l-e) , (A10) 
and, taking apocentre values and writing 

/-^, (All) 



we find 



l-e + 2 V " 1-e 
The two parts of the integral (A8) then give 



l+/ =k 1 + f) + T^~- (A12) 



L 



(A13) 



Since both x and \+ j (and 77) increase by 2-7T in one radial period we find that 
increases by 

L 

1 + 



$ = 7T 



< 2ir . (A14) 



Only when b = is $ = 2-7T and then both x and x+ reduce to 0. Thus the orbits 
do not close (unless y / L 2^ AGMb is rational) and they form rosettes with inner and 
outer radii of 



r p = ^Ja p (a p + 26) and r a = \J a a (a a + 2b) , (A15) 

where 

a a = a p (l + e)/(l - e) = (o - 6)(1 + e) . (A16) 

If in place of inventing \ an d X+ we proceed with (A4), we obtain the ugly 
formula (cf. (2.36K)) 



tan 1 I a I -tan n/2 | + ; = tan 1 J- - L tanr7/2| , 

1 W l + e /7 J VL 2 + 4GM6 VV 1 + e + / 



which is term for term the same as (A13). 
One may show from the above formulae that 



L 1-e 2 e. 



VL 2 + AGMb V 1 ~ e + e 



(A17) 



(A18) 



To get the general solution to the N-body problem one procedure is to take a 
value of 77, from that determine t from (A5), s from (A4) and r = V s 2 — b 2 and 
finally <p from (A17). We finally have 

r = r ( A sin <j) + B cos 0) 

as before. (A9) gives the basic relationship between r, s, 77, x and x+ whereas 
(A13) relates to % and x+ and (A5) relates t to 77. 
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Note Added in Proof 

Geometrical Insight gives greater clarity. From (2.5), the problem is a central 
one both in 3N dimensions and for each particle. Together the initial r and r 
define a plane in 3iV dimensions through the centre. The central force lies in 
that plane so all the motion continues in it. The motion is just planar motion 
under the potential V(r) so we get the usual orbit. The orbit of each particle is the 
projection of this 3N planar orbit onto the orbital plane of that particle (that has 
most of the 3 N coordinates fixed at zero) . Thus in all cases the particle orbits are 
projections of the 3iV planar orbit which is a familiar one. If this is an ellipse with 
the barycentre as its focus, then the particle orbit will be the projected ellipse 
but the projection does not preserve the barycentre as its focus. Similarly if the 
3iV orbit is a rosette, the individual particle orbits will be projected rosettes, i.e., 
rosettes between similar ellipses. That will happen for most potentials, e.g., the 
isochrone. 
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